#include "cellobject.h"

#include <cfloat>>

#define DTOR (pi/180.0f)

CellObject::CellObject()
{
    vertice_num = 0;
    translateX = translateY = translateZ = 0.0;
    rotateX = rotateZ = rotateY = 0.0;
    rotateCtrl = rotateAxisX = rotateAxisY = rotateAxisZ = 0.0;
    scaleX = scaleY = scaleZ = 1.0;
    is_surface_constructed = POINT_CLOUD;
    isCSG = false;
    AxisRotate();

    CellMath::identityMat(newRotateMatrix);
}

CellObject::~CellObject()
{
    //cout << "enter ~CellObject()" << endl;

    delete [] vertice; //*GLdouble
    delete [] normals; //*Gldouble
    delete [] face; //*int
    delete [] selectedFace; //*bool
    delete [] selectedPoint; //*bool
    delete [] deletePoint; //*bool
    delete [] deleteFace; //*bool
    if (intersectionPlane.size()>0) //vector<double>
        intersectionPlane.clear();

    //cout << "exit ~CellObject()" << endl;
}

CellObject* CellObject::clone()
{
    CellObject* obj = new CellObject();
    obj->vertice_num = this->vertice_num;
    obj->face_num = this->face_num;
    obj->type = this->type;
    obj->is_surface_constructed = this->is_surface_constructed;
    obj->vertice = new GLdouble[vertice_num*3];
    obj->normals = new GLdouble[vertice_num*3];
    obj->face = new int[face_num*3];
    obj->deletePoint = new bool[vertice_num];
    obj->selectedPoint = new bool[vertice_num];
    obj->deleteFace = new bool[face_num];
    obj->selectedFace = new bool[face_num];
    obj->translateX = this->translateX;
    obj->translateY = this->translateY;
    obj->translateZ = this->translateZ;
    obj->rotateX = this->rotateX;
    obj->rotateY = this->rotateY;
    obj->rotateZ = this->rotateZ;
    obj->rotateCtrl = this->rotateCtrl;
    obj->rotateAxisX = this->rotateAxisX;
    obj->rotateAxisY = this->rotateAxisY;
    obj->rotateAxisZ = this->rotateAxisZ;
    obj->scaleX = this->scaleX;
    obj->scaleY = this->scaleY;
    obj->scaleZ = this->scaleZ;
    obj->AxisRotate();

    for (int i=0;i<vertice_num*3;i++) {
        obj->vertice[i] = this->vertice[i];
        obj->normals[i] = this->normals[i];
    }
    for (int i=0;i<face_num*3;i++) {
        obj->face[i] = this->face[i];
    }
    for (int i=0;i<vertice_num;i++) {
        obj->deletePoint[i] = this->deletePoint[i];
        obj->selectedPoint[i] = this->selectedPoint[i];
    }
    for (int i=0;i<face_num;i++) {
        obj->deleteFace[i] = this->deleteFace[i];
        obj->selectedFace[i] = this->selectedFace[i];
    }
    obj->isCSG = this->isCSG;
    if (this->isCSG==true) {
        obj->csgObject = new CSG();
        obj->csgObject = this->csgObject->clone();
    }
    return obj;
}

void CellObject::AxisRotate()
{
    double matX[4][4], matY[4][4], matZ[4][4];

    //calculate the x, y, & z matrices based on the rotateN value
    CellMath::rotateMatX(rotateX, matX);
    CellMath::rotateMatY(rotateY, matY);
    CellMath::rotateMatZ(rotateZ, matZ);
    //calculate the rotation matrix by multiplying matZ*matY*matX in that order and store it in "result"
    double tmp[4][4], result[4][4];
    CellMath::matMutiply(matZ, matY, tmp);
    CellMath::matMutiply(tmp, matX, result);
    double tmpX[4], tmpY[4], tmpZ[4];
    /*
      tmpX = [1 0 0 1]
      tmpY = [0 1 0 1]
      tmpZ = [0 0 1 1]
    */
    tmpX[0] = tmpX[3] = 1.0;
    tmpX[1] = tmpX[2] = 0.0;
    tmpY[1] = tmpY[3] = 1.0;
    tmpY[0] = tmpY[2] = 0.0;
    tmpZ[2] = tmpZ[3] = 1.0;
    tmpZ[0] = tmpZ[1] = 0.0;
    CellMath::matMutiplyVector(result, tmpX, axisX);
    CellMath::matMutiplyVector(result, tmpY, axisY);
    CellMath::matMutiplyVector(result, tmpZ, axisZ);

    CellMath::rotateMatX(-rotateX, matX);
    CellMath::rotateMatY(-rotateY, matY);
    CellMath::rotateMatZ(-rotateZ, matZ);
    CellMath::matMutiply(matZ, matY, tmp);
    CellMath::matMutiply(tmp, matX, rotateMatrix);
}

void CellObject::printRotateMatrix()
{
    cout << rotateMatrix[0][0] << "\t" << rotateMatrix[0][1] << "\t" << rotateMatrix[0][2] << "\t" << rotateMatrix[0][3] << endl;
    cout << rotateMatrix[1][0] << "\t" << rotateMatrix[1][1] << "\t" << rotateMatrix[1][2] << "\t" << rotateMatrix[1][3] << endl;
    cout << rotateMatrix[2][0] << "\t" << rotateMatrix[2][1] << "\t" << rotateMatrix[2][2] << "\t" << rotateMatrix[2][3] << endl;
    cout << rotateMatrix[3][0] << "\t" << rotateMatrix[3][1] << "\t" << rotateMatrix[3][2] << "\t" << rotateMatrix[3][3] << endl << endl;
}

double CellObject::getMaxScale()
{
    if(scaleX >= scaleY && scaleX >= scaleZ)
        return scaleX;
    else if(scaleY >= scaleZ)
        return scaleY;
    else
        return scaleZ;
}

void CellObject::transformVertex()
{
    double matX[4][4], matY[4][4], matZ[4][4], matT[4][4], matS[4][4];
    CellMath::rotateMatX(rotateX, matX);
    CellMath::rotateMatY(rotateY, matY);
    CellMath::rotateMatZ(rotateZ, matZ);
    CellMath::translateMat(translateX, translateY, translateZ, matT);
    CellMath::scaleMat(scaleX, scaleY, scaleZ, matS);
    double tmp[4][4], result[4][4], matR[4][4], inverseResult[4][4], transpoteResult[4][4];
    CellMath::matMutiply(matZ, matY, tmp);
    CellMath::matMutiply(tmp, matX, matR);
    CellMath::matMutiply(matT, matR, tmp);
    CellMath::matMutiply(tmp, matS, result);
    CellMath::InverseMat(result, inverseResult);
    CellMath::TranspoteMat(inverseResult, transpoteResult);
    for (int i=0;i<vertice_num;i++) {
        double tmpV[4], tmpR[4];
        tmpV[0] = vertice[i*3+0];
        tmpV[1] = vertice[i*3+1];
        tmpV[2] = vertice[i*3+2];
        tmpV[3] = 1.0;
        CellMath::matMutiplyVector(result, tmpV, tmpR);
        vertice[i*3+0] = tmpR[0]/tmpR[3];
        vertice[i*3+1] = tmpR[1]/tmpR[3];
        vertice[i*3+2] = tmpR[2]/tmpR[3];
        tmpV[0] = normals[i*3+0];
        tmpV[1] = normals[i*3+1];
        tmpV[2] = normals[i*3+2];
        tmpV[3] = 1.0;
        CellMath::matMutiplyVector(transpoteResult, tmpV, tmpR);
        normals[i*3+0] = tmpR[0];
        normals[i*3+1] = tmpR[1];
        normals[i*3+2] = tmpR[2];
        double l = sqrt(normals[i*3+0]*normals[i*3+0]+normals[i*3+1]*normals[i*3+1]+normals[i*3+2]*normals[i*3+2]);
        normals[i*3+0] = normals[i*3+0]/l;
        normals[i*3+1] = normals[i*3+1]/l;
        normals[i*3+2] = normals[i*3+2]/l;
    }
}

void CellObject::transformCSG()
{
    double matX[4][4], matY[4][4], matZ[4][4], matT[4][4], matS[4][4];
    CellMath::rotateMatX(rotateX, matX);
    CellMath::rotateMatY(rotateY, matY);
    CellMath::rotateMatZ(rotateZ, matZ);
    CellMath::translateMat(translateX, translateY, translateZ, matT);
    CellMath::scaleMat(scaleX, scaleY, scaleZ, matS);
    double tmp[4][4], result[4][4], matR[4][4], inverseResult[4][4], transpoteResult[4][4];
    CellMath::matMutiply(matZ, matY, tmp);
    CellMath::matMutiply(tmp, matX, matR);
    CellMath::matMutiply(matT, matR, tmp);
    CellMath::matMutiply(tmp, matS, result);
    CellMath::InverseMat(result, inverseResult);
    CellMath::TranspoteMat(inverseResult, transpoteResult);
    for (size_t i=0;i<csgObject->polygons.size();i++) {
        for (size_t j=0;j<csgObject->polygons[i].vertices.size();j++) {
            double tmpV[4], tmpR[4];
            tmpV[0] = csgObject->polygons[i].vertices[j].pos.x;
            tmpV[1] = csgObject->polygons[i].vertices[j].pos.y;
            tmpV[2] = csgObject->polygons[i].vertices[j].pos.z;
            tmpV[3] = 1.0;
            CellMath::matMutiplyVector(result, tmpV, tmpR);
            csgObject->polygons[i].vertices[j].pos.x = tmpR[0]/tmpR[3];
            csgObject->polygons[i].vertices[j].pos.y = tmpR[1]/tmpR[3];
            csgObject->polygons[i].vertices[j].pos.z = tmpR[2]/tmpR[3];
            tmpV[0] = csgObject->polygons[i].vertices[j].normal.x;
            tmpV[1] = csgObject->polygons[i].vertices[j].normal.y;
            tmpV[2] = csgObject->polygons[i].vertices[j].normal.z;
            tmpV[3] = 1.0;
            CellMath::matMutiplyVector(transpoteResult, tmpV, tmpR);
            csgObject->polygons[i].vertices[j].normal.x = tmpR[0];
            csgObject->polygons[i].vertices[j].normal.y = tmpR[1];
            csgObject->polygons[i].vertices[j].normal.z = tmpR[2];
            csgObject->polygons[i].vertices[j].normal.Normalize();
        }
        csgObject->polygons[i].update();
    }
}

void CellObject::CreateCellObject(int num, vector<double> vec, vector<double> nor)
{
    type = CELL_MESH;
    if (vertice_num!=0) {
        delete [] vertice;
        delete [] normals;
        delete [] selectedPoint;
    }
    vertice = new double[num*3];
    normals = new double[num*3];
    selectedPoint = new bool[num];
    deletePoint = new bool[num];
    for (int i=0;i<num;i++)
        deletePoint[i]=false;
    for (int i=0;i<num;i++) {
        vertice[i*3+0] = vec[i*3+0]/100.0;
        vertice[i*3+1] = vec[i*3+1]/100.0;
        vertice[i*3+2] = vec[i*3+2]/100.0;
        normals[i*3+0] = nor[i*3+0];
        normals[i*3+1] = nor[i*3+1];
        normals[i*3+2] = nor[i*3+2];
    }
    vertice_num = num;
    resetSelectedPoint();
    //SurfaceReconstruction();
}

void CellObject::SaveCellObject(char *filename)
{
    FILE *out;
    out = fopen(filename, "w+");
    fprintf(out, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n");
    fprintf(out, "<MultiCellXML>\n");
    fprintf(out, "\t<BasementMembrane>\n");
    fprintf(out, "\t\t<BasementMembraneNodes>\n");
    for (int i=0;i<vertice_num;i++) {
        fprintf(out, "\t\t\t<BasementMembraneNode ID=\"%d\">\n", i);
        fprintf(out, "\t\t\t\t<Position units=\"microns\">(%lf,%lf,%lf)</Position>\n",
                vertice[3*i+0], vertice[3*i+1], vertice[3*i+2]);
        fprintf(out, "\t\t\t\t<Normal units=\"microns\">(%lf,%lf,%lf)</Normal>\n",
                normals[3*i+0], normals[3*i+1], normals[3*i+2]);
        fprintf(out, "\t\t\t</BasementMembraneNode>\n");
    }
    fprintf(out, "\t\t</BasementMembraneNodes>\n");


    fprintf(out, "\t\t<BasementMembraneFacets>\n");
    for (int i=0;i<face_num;i++) {
        fprintf(out, "\t\t\t<BasementMembraneFacet ID=\"%d\">\n", i);
        fprintf(out, "\t\t\t\t<Index>(%d,%d,%d)</Index>\n",
                face[3*i+0], face[3*i+1], face[3*i+2]);
        fprintf(out, "\t\t\t</BasementMembraneFacet>\n");
    }
    fprintf(out, "\t\t</BasementMembraneFacets>\n");
    fprintf(out, "\t</BasementMembrane>\n");
    fprintf(out, "</MultiCellXML>");
    fclose(out);
}

void CellObject::DrawCellObject(int id, int selected, int displayMode)
{
    if (vertice_num>0)
    {
        glPushMatrix();
        glLoadName(id);
        glTranslated(translateX, translateY, translateZ);
        glRotated(rotateCtrl, rotateAxisX, rotateAxisY, rotateAxisZ);
//        glRotated(rotateZ, 0.0, 0.0, 1.0);
//        glRotated(rotateY, 0.0, 1.0, 0.0);
//        glRotated(rotateX, 1.0, 0.0, 0.0);
        glScaled(scaleX, scaleY, scaleZ);

        for (int i=0; i<vertice_num*3; i+=3)
        {
            double vertex[4] = {vertice[i], vertice[i+1], vertice[i+2], 1};
            double result[4];

            CellMath::matMutiplyVector(newRotateMatrix, vertex, result);
            transformedVertices[i] = result[0];
            transformedVertices[i+1] = result[1];
            transformedVertices[i+2] = result[2];
        }

        if (is_surface_constructed == POINT_CLOUD)
        {
            if (selected!=id)
                glColor3f(0.0,0.0,1.0);
            else
                glColor3f(1.0,1.0,0.0);
            glPointSize(8.0);
            glEnable(GL_POINT_SMOOTH);
            glEnable(GL_BLEND);
            glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
            glBegin(GL_POINTS);
            for (int i=0;i<vertice_num;i++) {
                glVertex3d(transformedVertices[i*3], transformedVertices[i*3+1], transformedVertices[i*3+2]);
            }

            glEnd();
        }
        else if (is_surface_constructed == TRIANGLES)
        {
            if (displayMode==CELL_TRIANGLE)
                glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
            else if (displayMode==CELL_LINE)
                glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
            else {
                glPolygonMode(GL_FRONT_AND_BACK, GL_POINT);
                glPointSize(2.0);
            }
            if (selected!=id) {
                GLfloat dif[4] = { 0.2, 0.2, 1.0, 1.0 };
                GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
                GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
                glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
                glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
                glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
                glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);
            } else {
                GLfloat dif[4] = { 0.2, 1.0, 0.2, 1.0 };
                GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
                GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
                glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
                glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
                glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
                glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);
            }
            glEnable(GL_LIGHTING);
            glDrawElements(GL_TRIANGLES, face_num*3, GL_UNSIGNED_INT, face);
            glDisable(GL_LIGHTING);
            glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
        }
        else if (is_surface_constructed == LINE_SEGMENT)
        {
            if (selected!=id)
            {
                glColor3f(0.2,0.2,1.0);
            }
            else
            {
                glColor3f(0.2,1.0,0.2);
            }
            glEnable(GL_LINE_SMOOTH);
            glLineWidth(5.0);
            glBegin(GL_LINES);
            glVertex3d(vertice[0], vertice[1], vertice[2]);
            glVertex3d(vertice[3], vertice[4], vertice[5]);
            glEnd();
        }
        glPopMatrix();
    }
}

void CellObject::DrawCellObjectDirect(int id, int selected, int secondSelected, int displayMode)
{
    if (vertice_num>0) {
        glPushMatrix();
        glLoadName(id);
        glTranslated(translateX, translateY, translateZ);
        glRotated(rotateCtrl, rotateAxisX, rotateAxisY, rotateAxisZ);
//        glRotated(rotateZ, 0.0, 0.0, 1.0);
//        glRotated(rotateY, 0.0, 1.0, 0.0);
//        glRotated(rotateX, 1.0, 0.0, 0.0);
        glScaled(scaleX, scaleY, scaleZ);

        glEnable(GL_BLEND);
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

        if (is_surface_constructed == POINT_CLOUD)
        {
            if (selected!=id)
                glColor3f(0.2,0.2,1.0);
            else
                glColor3f(0.2,1.0,0.2);
            glPointSize(8.0);
            glEnable(GL_POINT_SMOOTH);
            glBegin(GL_POINTS);
            for (int i=0;i<vertice_num;i++) {
                if (deletePoint[i]==false) {
                    glVertex3d(vertice[i*3], vertice[i*3+1], vertice[i*3+2]);
                }
            }
            glEnd();
        }
        else if(is_surface_constructed == TRIANGLES)
        {
            if (displayMode==CELL_TRIANGLE)
                glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
            else if (displayMode==CELL_LINE)
                glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
            else {
                glPolygonMode(GL_FRONT_AND_BACK, GL_POINT);
                glPointSize(4.0);
                glEnable(GL_POINT_SMOOTH);
            }
            if (selected==id)
            {
                GLfloat dif[4] = { 0.2, 1.0, 0.2, 1.0 };
                GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
                GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
                glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
                glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
                glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
                glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);
            }
            else if (secondSelected==id)
            {
                GLfloat dif[4] = { 1.0, 0.2, 0.2, 1.0 };
                GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
                GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
                glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
                glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
                glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
                glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);
            }
            else
            {
                //GLfloat dif[4] = { 0.2, 0.2, 1.0, 1.0 };
                GLfloat dif[4] = { 0.2, 0.2, 1.0, 0.7 };
                GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
                GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
                glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
                glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
                glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
                glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);
            }
            glEnable(GL_LIGHTING);

            //apply rotation matrix to the set of vertice and store it in another array in order to keep a copy of the original shape
            for (int i=0; i<vertice_num*3; i+=3)
            {
                double vertex[4] = {vertice[i], vertice[i+1], vertice[i+2], 1};
                double result[4];

                CellMath::matMutiplyVector(newRotateMatrix, vertex, result);
                transformedVertices[i] = result[0];
                transformedVertices[i+1] = result[1];
                transformedVertices[i+2] = result[2];
            }

            //draw one triangle at a time using the 3d coordinates and normals
            for (int i=0;i<face_num;i++)
            {
                if (deleteFace[i]==false)
                {
                    glBegin(GL_TRIANGLES);
                    glNormal3d(normals[face[i*3+0]*3+0], normals[face[i*3+0]*3+1], normals[face[i*3+0]*3+2]);
                    glVertex3d(vertice[face[i*3+0]*3+0], vertice[face[i*3+0]*3+1], vertice[face[i*3+0]*3+2]);
                    //glVertex3d(transformedVertices[face[i*3+0]*3+0], transformedVertices[face[i*3+0]*3+1], transformedVertices[face[i*3+0]*3+2]);
                    glNormal3d(normals[face[i*3+1]*3+0], normals[face[i*3+1]*3+1], normals[face[i*3+1]*3+2]);
                    glVertex3d(vertice[face[i*3+1]*3+0], vertice[face[i*3+1]*3+1], vertice[face[i*3+1]*3+2]);
                    //glVertex3d(transformedVertices[face[i*3+1]*3+0], transformedVertices[face[i*3+1]*3+1], transformedVertices[face[i*3+1]*3+2]);
                    glNormal3d(normals[face[i*3+2]*3+0], normals[face[i*3+2]*3+1], normals[face[i*3+2]*3+2]);
                    glVertex3d(vertice[face[i*3+2]*3+0], vertice[face[i*3+2]*3+1], vertice[face[i*3+2]*3+2]);
                    //glVertex3d(transformedVertices[face[i*3+2]*3+0], transformedVertices[face[i*3+2]*3+1], transformedVertices[face[i*3+2]*3+2]);
                    glEnd();
                }
            }
            glDisable(GL_LIGHTING);
            glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
        }
        else if (is_surface_constructed == LINE_SEGMENT)
        {
            if (selected!=id)
            {
                glColor3f(0.2,0.2,1.0);
            }
            else
            {
                glColor3f(0.2,1.0,0.2);
            }
            glEnable(GL_LINE_SMOOTH);
            glLineWidth(5.0);
            glBegin(GL_LINES);
            glVertex3d(vertice[0], vertice[1], vertice[2]);
            glVertex3d(vertice[3], vertice[4], vertice[5]);
            glEnd();
        }
        glPopMatrix();
    }
}

void CellObject::DrawCellObjectDirectTransparent(int id, int displayMode, int transparencyLvl, int color)
{
    if (vertice_num>0) {
        //setup
        glPushMatrix();
        glLoadName(id);
        glTranslated(translateX, translateY, translateZ);
        glRotated(rotateCtrl, rotateAxisX, rotateAxisY, rotateAxisZ);
//        glRotated(rotateZ, 0.0, 0.0, 1.0);
//        glRotated(rotateY, 0.0, 1.0, 0.0);
//        glRotated(rotateX, 1.0, 0.0, 0.0);
        glScaled(scaleX, scaleY, scaleZ);

        //triangle mode
        glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);

        glEnable(GL_BLEND);
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

        //PROCESS PARAMETERS
        double transparent;
        if(transparencyLvl == LOW_TRANSPARENCY)
            transparent = 0.6;
        else if(transparencyLvl == MID_TRANSPARENCY)
            transparent = 0.4;
        else //HIGH_TRANSPARENCY
            transparent = 0.2;

        double red, green, blue;
        if(color == PURPLE) {
            red = 1.0;
            green = 0.2;
            blue = 1.0;
        } else {
            red = 1.0;
            green = 1.0;
            blue = 1.0;
        }

        //GL_Diffuse = { R, G, B, Transparency}
        GLfloat dif[4] = { red, green, blue, transparent };
        GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
        GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
        glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
        glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
        glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
        glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);

        glEnable(GL_LIGHTING);

        for (int i=0; i<vertice_num*3; i+=3)
        {
            double vertex[4] = {vertice[i], vertice[i+1], vertice[i+2], 1};
            double result[4];

            CellMath::matMutiplyVector(newRotateMatrix, vertex, result);
            transformedVertices[i] = result[0];
            transformedVertices[i+1] = result[1];
            transformedVertices[i+2] = result[2];
        }

        if (is_surface_constructed == POINT_CLOUD)
        {
            glPointSize(8.0);
            glEnable(GL_POINT_SMOOTH);
            glEnable(GL_BLEND);
            glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
            glBegin(GL_POINTS);
            for (int i=0;i<vertice_num;i++) {
                glVertex3d(transformedVertices[i*3], transformedVertices[i*3+1], transformedVertices[i*3+2]);
            }

            glEnd();
        }
        else if (is_surface_constructed == TRIANGLES)
        {
            for (int i=0;i<face_num;i++)
            {
                if (deleteFace[i]==false)
                {
                    glBegin(GL_TRIANGLES);
                    glNormal3d(normals[face[i*3+0]*3+0], normals[face[i*3+0]*3+1], normals[face[i*3+0]*3+2]);
                    glVertex3d(transformedVertices[face[i*3+0]*3+0], transformedVertices[face[i*3+0]*3+1], transformedVertices[face[i*3+0]*3+2]);
                    glNormal3d(normals[face[i*3+1]*3+0], normals[face[i*3+1]*3+1], normals[face[i*3+1]*3+2]);
                    glVertex3d(transformedVertices[face[i*3+1]*3+0], transformedVertices[face[i*3+1]*3+1], transformedVertices[face[i*3+1]*3+2]);
                    glNormal3d(normals[face[i*3+2]*3+0], normals[face[i*3+2]*3+1], normals[face[i*3+2]*3+2]);
                    glVertex3d(transformedVertices[face[i*3+2]*3+0], transformedVertices[face[i*3+2]*3+1], transformedVertices[face[i*3+2]*3+2]);
                    glEnd();
                }
            }
        }
        else if (is_surface_constructed == LINE_SEGMENT)
        {
            glEnable(GL_LINE_SMOOTH);
            glLineWidth(5.0);
            glBegin(GL_LINES);
            glVertex3d(vertice[0], vertice[1], vertice[2]);
            glVertex3d(vertice[3], vertice[4], vertice[5]);
            glEnd();
        }

        glDisable(GL_LIGHTING);
        glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);

        glPopMatrix();
    }
}

void CellObject::SurfaceReconstruction(double radius)
{
    BPA *bpa = new BPA(vertice, normals, vertice_num, radius);
    bpa->Reconstruction(face_num);
    face = new int[face_num*3];
    selectedFace = new bool[face_num];
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    resetSelectedFace();
    bpa->CopyFaceIndex(face);
    is_surface_constructed = TRIANGLES;
    glEnableClientState(GL_VERTEX_ARRAY);
    glEnableClientState(GL_NORMAL_ARRAY);
    glNormalPointer(GL_DOUBLE, 0, normals);
    glVertexPointer(3, GL_DOUBLE, 0, vertice);
}

bool CellObject::isSurfaceConstructed()
{
    if (is_surface_constructed == TRIANGLES)
    {
        return true;
    }
    else
    {
        return false;
    }
}

void CellObject::ComputeIntersectionPlane(double a, double b, double c, double d)
{
    if (is_surface_constructed == TRIANGLES) {
        intersectionPlane.clear();
        double matX[4][4], matY[4][4], matZ[4][4], matT[4][4], matS[4][4];
        CellMath::rotateMatX(rotateX, matX);
        CellMath::rotateMatY(rotateY, matY);
        CellMath::rotateMatZ(rotateZ, matZ);
        CellMath::translateMat(translateX, translateY, translateZ, matT);
        CellMath::scaleMat(scaleX, scaleY, scaleZ, matS);
        double tmp[4][4], result[4][4], matR[4][4], inverseResult[4][4], transpoteResult[4][4];
        CellMath::matMutiply(matZ, matY, tmp);
        CellMath::matMutiply(tmp, matX, matR);
        CellMath::matMutiply(matT, matR, tmp);
        CellMath::matMutiply(tmp, matS, result);
        CellMath::InverseMat(result, inverseResult);
        CellMath::TranspoteMat(inverseResult, transpoteResult);
        for (int i=0;i<face_num;i++) {
            double tmpV0[4], tmpR0[4], tmpV1[4], tmpR1[4], tmpV2[4], tmpR2[4];
            tmpV0[0] = vertice[face[i*3+0]*3+0];
            tmpV0[1] = vertice[face[i*3+0]*3+1];
            tmpV0[2] = vertice[face[i*3+0]*3+2];
            tmpV0[3] = 1.0;
            CellMath::matMutiplyVector(result, tmpV0, tmpR0);
            tmpR0[0]/=tmpR0[3];
            tmpR0[1]/=tmpR0[3];
            tmpR0[2]/=tmpR0[3];
            tmpV1[0] = vertice[face[i*3+1]*3+0];
            tmpV1[1] = vertice[face[i*3+1]*3+1];
            tmpV1[2] = vertice[face[i*3+1]*3+2];
            tmpV1[3] = 1.0;
            CellMath::matMutiplyVector(result, tmpV1, tmpR1);
            tmpR1[0]/=tmpR1[3];
            tmpR1[1]/=tmpR1[3];
            tmpR1[2]/=tmpR1[3];
            tmpV2[0] = vertice[face[i*3+2]*3+0];
            tmpV2[1] = vertice[face[i*3+2]*3+1];
            tmpV2[2] = vertice[face[i*3+2]*3+2];
            tmpV2[3] = 1.0;
            CellMath::matMutiplyVector(result, tmpV2, tmpR2);
            tmpR2[0]/=tmpR2[3];
            tmpR2[1]/=tmpR2[3];
            tmpR2[2]/=tmpR2[3];
            double rvx0, rvy0, rvz0, rvx1, rvy1, rvz1;
            bool flag = CellMath::triangleIntersection(a, b, c, d,
                                                       tmpR0[0], tmpR0[1], tmpR0[2],
                    tmpR1[0], tmpR1[1], tmpR1[2],
                    tmpR2[0], tmpR2[1], tmpR2[2],
                    rvx0, rvy0, rvz0, rvx1, rvy1, rvz1);
            if (flag) {
                intersectionPlane.push_back(rvx0);
                intersectionPlane.push_back(rvy0);
                intersectionPlane.push_back(rvz0);
                intersectionPlane.push_back(rvx1);
                intersectionPlane.push_back(rvy1);
                intersectionPlane.push_back(rvz1);
            }
        }
    }
}

void CellObject::DrawIntersectionPlane()
{
    glPushMatrix();
    glEnable(GL_LINE_SMOOTH);
    glEnable(GL_BLEND);
    glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
    glHint(GL_LINE_SMOOTH_HINT, GL_DONT_CARE);
    glLineWidth(4.0);
    glBegin(GL_LINES);
    for (size_t i=0; i<intersectionPlane.size()/6; i++) {
        glVertex3d(intersectionPlane[i*6+0], intersectionPlane[i*6+1], intersectionPlane[i*6+2]);
        glVertex3d(intersectionPlane[i*6+3], intersectionPlane[i*6+4], intersectionPlane[i*6+5]);
    }
    glEnd();
    glPopMatrix();
}

void CellObject::drawEditPoint()
{
    glPushMatrix();
    glPointSize(2.0);
    glTranslated(translateX, translateY, translateZ);
    glRotated(rotateZ, 0.0, 0.0, 1.0);
    glRotated(rotateY, 0.0, 1.0, 0.0);
    glRotated(rotateX, 1.0, 0.0, 0.0);
    glScaled(scaleX, scaleY, scaleZ);
    for (int i=0;i<vertice_num;i++) {
        glPushMatrix();
        glLoadName(i);
        glPointSize(8.0);
        glEnable(GL_POINT_SMOOTH);
        glEnable(GL_BLEND);
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
        glBegin(GL_POINTS);
        if (selectedPoint[i]==true) {
            glColor3f(1.0,0.2,0.2);
            if (deletePoint[i]==false) {
                glVertex3d(vertice[i*3], vertice[i*3+1], vertice[i*3+2]);
            }
        } else {
            glColor3f(0.5,0.5,0.5);
            if (deletePoint[i]==false) {
                glVertex3d(vertice[i*3], vertice[i*3+1], vertice[i*3+2]);
            }
        }
        glEnd();
        glPopMatrix();
    }
    glPopMatrix();
}

void CellObject::drawEditFaces()
{
    glPushMatrix();
    glTranslated(translateX, translateY, translateZ);
    glRotated(rotateZ, 0.0, 0.0, 1.0);
    glRotated(rotateY, 0.0, 1.0, 0.0);
    glRotated(rotateX, 1.0, 0.0, 0.0);
    glScaled(scaleX, scaleY, scaleZ);
    glEnable(GL_LIGHTING);
    for (int i=0;i<face_num;i++) {
        glPushMatrix();
        glLoadName(i);
        glBegin(GL_TRIANGLES);
        if (selectedFace[i]==true) {
            GLfloat dif[4] = { 1.0, 0.2, 0.2, 1.0 };
            GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
            GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
            glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
            glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
            glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
            glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);
            if (deleteFace[i]==false) {
                glNormal3d(normals[face[i*3+0]*3+0], normals[face[i*3+0]*3+1], normals[face[i*3+0]*3+2]);
                glVertex3d(vertice[face[i*3+0]*3+0], vertice[face[i*3+0]*3+1], vertice[face[i*3+0]*3+2]);
                glNormal3d(normals[face[i*3+1]*3+0], normals[face[i*3+1]*3+1], normals[face[i*3+1]*3+2]);
                glVertex3d(vertice[face[i*3+1]*3+0], vertice[face[i*3+1]*3+1], vertice[face[i*3+1]*3+2]);
                glNormal3d(normals[face[i*3+2]*3+0], normals[face[i*3+2]*3+1], normals[face[i*3+2]*3+2]);
                glVertex3d(vertice[face[i*3+2]*3+0], vertice[face[i*3+2]*3+1], vertice[face[i*3+2]*3+2]);
            }
        } else {
            GLfloat dif[4] = { 0.5, 0.5, 0.5, 1.0 };
            GLfloat amb[4] = { 0.0, 0.0, 0.0, 1.0 };
            GLfloat spec[4] = { 0.774, 0.774, 0.774, 1.0 };
            glMaterialfv(GL_FRONT, GL_DIFFUSE, dif);
            glMaterialfv(GL_FRONT, GL_AMBIENT, amb);
            glMaterialfv(GL_FRONT, GL_SPECULAR, spec);
            glMaterialf(GL_FRONT, GL_SHININESS, 0.6*128.0);
            if (deleteFace[i]==false) {
                glNormal3d(normals[face[i*3+0]*3+0], normals[face[i*3+0]*3+1], normals[face[i*3+0]*3+2]);
                glVertex3d(vertice[face[i*3+0]*3+0], vertice[face[i*3+0]*3+1], vertice[face[i*3+0]*3+2]);
                glNormal3d(normals[face[i*3+1]*3+0], normals[face[i*3+1]*3+1], normals[face[i*3+1]*3+2]);
                glVertex3d(vertice[face[i*3+1]*3+0], vertice[face[i*3+1]*3+1], vertice[face[i*3+1]*3+2]);
                glNormal3d(normals[face[i*3+2]*3+0], normals[face[i*3+2]*3+1], normals[face[i*3+2]*3+2]);
                glVertex3d(vertice[face[i*3+2]*3+0], vertice[face[i*3+2]*3+1], vertice[face[i*3+2]*3+2]);
            }
        }
        glEnd();
        glPopMatrix();
    }
    glDisable(GL_LIGHTING);
    glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
    glPopMatrix();
}

void CellObject::resetSelectedPoint()
{
    for (int i=0;i<vertice_num;i++)
        selectedPoint[i]=false;
}

void CellObject::setPointSelect(int index){
    if (index<vertice_num&&index>=0)
        selectedPoint[index]=true;
}

void CellObject::resetSelectedFace()
{
    for (int i=0;i<face_num;i++)
        selectedFace[i]=false;
}

void CellObject::setFaceSelect(int index){
    if (index<face_num&&index>=0)
        selectedFace[index]=true;
}

void CellObject::setPointsSelectInArea(Vector3 point[8], Vector3 view, Vector3 up)
{
    for (int i=0; i<vertice_num; i++) {
        Vector3 p = Vector3(vertice[i*3], vertice[i*3+1], vertice[i*3+2]);
        if (CellMath::isVectorInFrustum(p, point, view, up))
            setPointSelect(i);
    }
}

void CellObject::setFacesSelectInArea(Vector3 point[8], Vector3 view, Vector3 up)
{
    for (int i=0; i<face_num; i++) {
        Vector3 p0 = Vector3(vertice[face[i*3+0]*3+0], vertice[face[i*3+0]*3+1], vertice[face[i*3+0]*3+2]);
        Vector3 p1 = Vector3(vertice[face[i*3+1]*3+0], vertice[face[i*3+1]*3+1], vertice[face[i*3+1]*3+2]);
        Vector3 p2 = Vector3(vertice[face[i*3+2]*3+0], vertice[face[i*3+2]*3+1], vertice[face[i*3+2]*3+2]);
        if (CellMath::isVectorInFrustum(p0, point, view, up)&&CellMath::isVectorInFrustum(p1, point, view, up)&&
                CellMath::isVectorInFrustum(p2, point, view, up))
            setFaceSelect(i);
    }
}

int CellObject::selectPointsCenter(Vector3 &center)
{
    int num=0;
    center = Vector3(0.0, 0.0, 0.0);
    for (int i=0;i<vertice_num;i++) {
        if (selectedPoint[i]==true) {
            num++;
            center = center+Vector3(vertice[i*3], vertice[i*3+1], vertice[i*3+2]);
        }
    }
    if (num>0) {
        center=center/num;
        double matX[4][4], matY[4][4], matZ[4][4], matT[4][4], matS[4][4];
        CellMath::rotateMatX(rotateX, matX);
        CellMath::rotateMatY(rotateY, matY);
        CellMath::rotateMatZ(rotateZ, matZ);
        CellMath::translateMat(translateX, translateY, translateZ, matT);
        CellMath::scaleMat(scaleX, scaleY, scaleZ, matS);
        double tmp[4][4], result[4][4], matR[4][4];
        CellMath::matMutiply(matZ, matY, tmp);
        CellMath::matMutiply(tmp, matX, matR);
        CellMath::matMutiply(matT, matR, tmp);
        CellMath::matMutiply(tmp, matS, result);

        double tmpV[4], tmpR[4];
        tmpV[0] = center.x;
        tmpV[1] = center.y;
        tmpV[2] = center.z;
        tmpV[3] = 1.0;
        CellMath::matMutiplyVector(result, tmpV, tmpR);
        center.x = tmpR[0]/tmpR[3];
        center.y = tmpR[1]/tmpR[3];
        center.z = tmpR[2]/tmpR[3];
    }
    return num;
}

void CellObject::moveSelectPoints(Vector3 vec)
{
    for (int i=0;i<vertice_num;i++) {
        if (selectedPoint[i]==true) {
            vertice[i*3]+=vec.x;
            vertice[i*3+1]+=vec.y;
            vertice[i*3+2]+=vec.z;
        }
    }
}

void CellObject::moveSelectPointsFeather(Vector3 vec, double percent)
{
    int num=0;
    Vector3 center = Vector3(0.0, 0.0, 0.0);
    Vector3 aabb_min = Vector3(10000000.0, 10000000.0, 10000000.0);
    Vector3 aabb_max = Vector3(-10000000.0, -10000000.0, -10000000.0);
    for (int i=0;i<vertice_num;i++) {
        if (selectedPoint[i]==true) {
            num++;
            center = center+Vector3(vertice[i*3], vertice[i*3+1], vertice[i*3+2]);
            aabb_min = Vector3(min(aabb_min.x, vertice[i*3]), min(aabb_min.y, vertice[i*3+1]), min(aabb_min.z, vertice[i*3+2]));
            aabb_max = Vector3(max(aabb_max.x, vertice[i*3]), max(aabb_max.y, vertice[i*3+1]), max(aabb_max.z, vertice[i*3+2]));
        }
    }
    if (num>0) {
        if (num==1) {
            for (int i=0;i<vertice_num;i++) {
                if (selectedPoint[i]==true) {
                    vertice[i*3]+=vec.x;
                    vertice[i*3+1]+=vec.y;
                    vertice[i*3+2]+=vec.z;
                }
            }
        } else {
            center=center/num;
            double radius = (center-aabb_min).Lenght();
            for (int i=0;i<vertice_num;i++) {
                if (selectedPoint[i]==true) {
                    vertice[i*3]+=vec.x;
                    vertice[i*3+1]+=vec.y;
                    vertice[i*3+2]+=vec.z;
                } else {
                    double l = (Vector3(vertice[i*3], vertice[i*3+1], vertice[i*3+2])-center).Lenght();
                    if (l<radius*(percent+1.0)) {
                        vertice[i*3]+=vec.x*min(1.0, (1.0-(l-radius)/(radius*percent)));
                        vertice[i*3+1]+=vec.y*min(1.0, (1.0-(l-radius)/(radius*percent)));
                        vertice[i*3+2]+=vec.z*min(1.0, (1.0-(l-radius)/(radius*percent)));
                    }
                }
            }
        }
    }
}

int CellObject::selectFacesCenter(Vector3 &center)
{
    int num=0;
    center = Vector3(0.0, 0.0, 0.0);
    for (int i=0;i<face_num;i++) {
        if (selectedFace[i]==true) {
            num++;
            center = center+Vector3(vertice[face[i*3+0]*3+0], vertice[face[i*3+0]*3+1], vertice[face[i*3+0]*3+2])+
                    Vector3(vertice[face[i*3+1]*3+0], vertice[face[i*3+1]*3+1], vertice[face[i*3+1]*3+2])+
                    Vector3(vertice[face[i*3+2]*3+0], vertice[face[i*3+2]*3+1], vertice[face[i*3+2]*3+2]);
        }
    }
    if (num>0) {
        center=center/(num*3);
        double matX[4][4], matY[4][4], matZ[4][4], matT[4][4], matS[4][4];
        CellMath::rotateMatX(rotateX, matX);
        CellMath::rotateMatY(rotateY, matY);
        CellMath::rotateMatZ(rotateZ, matZ);
        CellMath::translateMat(translateX, translateY, translateZ, matT);
        CellMath::scaleMat(scaleX, scaleY, scaleZ, matS);
        double tmp[4][4], result[4][4], matR[4][4];
        CellMath::matMutiply(matZ, matY, tmp);
        CellMath::matMutiply(tmp, matX, matR);
        CellMath::matMutiply(matT, matR, tmp);
        CellMath::matMutiply(tmp, matS, result);

        double tmpV[4], tmpR[4];
        tmpV[0] = center.x;
        tmpV[1] = center.y;
        tmpV[2] = center.z;
        tmpV[3] = 1.0;
        CellMath::matMutiplyVector(result, tmpV, tmpR);
        center.x = tmpR[0]/tmpR[3];
        center.y = tmpR[1]/tmpR[3];
        center.z = tmpR[2]/tmpR[3];
    }
    return num;
}

void CellObject::moveSelectFaces(Vector3 vec)
{
    resetSelectedPoint();
    for (int i=0;i<face_num;i++) {
        if (selectedFace[i]==true) {
            setPointSelect(face[i*3+0]);
            setPointSelect(face[i*3+1]);
            setPointSelect(face[i*3+2]);
        }
    }
    for (int i=0;i<vertice_num;i++) {
        if (selectedPoint[i]==true) {
            vertice[i*3]+=vec.x;
            vertice[i*3+1]+=vec.y;
            vertice[i*3+2]+=vec.z;
        }
    }
    resetSelectedPoint();
}

void CellObject::moveSelectFacesFeather(Vector3 vec, double percent)
{
    resetSelectedPoint();
    for (int i=0;i<face_num;i++) {
        if (selectedFace[i]==true) {
            setPointSelect(face[i*3+0]);
            setPointSelect(face[i*3+1]);
            setPointSelect(face[i*3+2]);
        }
    }
    int num=0;
    Vector3 center = Vector3(0.0, 0.0, 0.0);
    Vector3 aabb_min = Vector3(10000000.0, 10000000.0, 10000000.0);
    Vector3 aabb_max = Vector3(-10000000.0, -10000000.0, -10000000.0);
    for (int i=0;i<vertice_num;i++) {
        if (selectedPoint[i]==true) {
            num++;
            center = center+Vector3(vertice[i*3], vertice[i*3+1], vertice[i*3+2]);
            aabb_min = Vector3(min(aabb_min.x, vertice[i*3]), min(aabb_min.y, vertice[i*3+1]), min(aabb_min.z, vertice[i*3+2]));
            aabb_max = Vector3(max(aabb_max.x, vertice[i*3]), max(aabb_max.y, vertice[i*3+1]), max(aabb_max.z, vertice[i*3+2]));
        }
    }
    if (num>0) {
        if (num==1) {
            for (int i=0;i<vertice_num;i++) {
                if (selectedPoint[i]==true) {
                    vertice[i*3]+=vec.x;
                    vertice[i*3+1]+=vec.y;
                    vertice[i*3+2]+=vec.z;
                }
            }
        } else {
            center=center/num;
            double radius = (center-aabb_min).Lenght();
            for (int i=0;i<vertice_num;i++) {
                if (selectedPoint[i]==true) {
                    vertice[i*3]+=vec.x;
                    vertice[i*3+1]+=vec.y;
                    vertice[i*3+2]+=vec.z;
                } else {
                    double l = (Vector3(vertice[i*3], vertice[i*3+1], vertice[i*3+2])-center).Lenght();
                    if (l<radius*(percent+1.0)) {
                        vertice[i*3]+=vec.x*min(1.0, (1.0-(l-radius)/(radius*percent)));
                        vertice[i*3+1]+=vec.y*min(1.0, (1.0-(l-radius)/(radius*percent)));
                        vertice[i*3+2]+=vec.z*min(1.0, (1.0-(l-radius)/(radius*percent)));
                    }
                }
            }
        }
    }
    resetSelectedPoint();
}

void CellObject::deleteSelectPoint()
{
    for (int i=0;i<vertice_num;i++) {
        if (selectedPoint[i]==true)
            deletePoint[i]=true;
    }
    for (int i=0;i<face_num;i++) {
        if (selectedPoint[face[i*3+0]]==true||selectedPoint[face[i*3+1]]==true||selectedPoint[face[i*3+2]]==true)
            deleteFace[i]=true;
    }
}

void CellObject::deleteSelectFace()
{
    bool *remainPoint = new bool[vertice_num];
    for (int i=0;i<vertice_num;i++)
        remainPoint[i] = false;
    for (int i=0;i<face_num;i++) {
        if (selectedFace[i]==true)
            deleteFace[i]=true;
        else {
            remainPoint[face[i*3+0]]=true;
            remainPoint[face[i*3+1]]=true;
            remainPoint[face[i*3+2]]=true;
        }
    }
    for (int i=0;i<vertice_num;i++) {
        if (remainPoint[i]==false)
            deletePoint[i]=true;
    }
    delete [] remainPoint;
}

void CellObject::buildSphere(double radius, int hseg, int wseg)
{
    type =CELL_SPHERE;
    is_surface_constructed = TRIANGLES;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }
    double dphi, dtheta;
    if (radius==0.0)
        radius = 5.0;
    sphereRadius = radius;
    if (hseg==0) {
        dphi = 15;
        sphereHeightSeg = 12;
    }
    else {
        dphi = 180.0/hseg;
        sphereHeightSeg = hseg;
    }
    if (wseg==0) {
        dtheta = 15;
        sphereAxisSeg = 24;
    }
    else {
        dtheta = 360.0/wseg;
        sphereAxisSeg = wseg;
    }

    vertice_num = sphereAxisSeg*(sphereHeightSeg-1)+2;
    face_num = sphereAxisSeg*(sphereHeightSeg-2)*2+sphereAxisSeg*2;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];
    double theta, phi;

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    int n = 0;
    vertice[3*n+0] = 0;
    vertice[3*n+1] = radius;
    vertice[3*n+2] = 0;
    normals[3*n+0] = 0;
    normals[3*n+1] = 1.0;
    normals[3*n+2] = 0;
    n++;
    for (phi=90-dphi; phi>-90.0+CELL_EPS; phi-=dphi)
        for (theta=0; theta<360-CELL_EPS; theta += dtheta) {
            vertice[3*n+0] = radius * cosf(phi*DTOR) * cosf(DTOR*theta);
            vertice[3*n+1] = radius * sinf(phi*DTOR);
            vertice[3*n+2] = radius * cosf(phi*DTOR) * sinf(DTOR*theta);
            Vector3 nor = Vector3(vertice[3*n+0], vertice[3*n+1], vertice[3*n+2]);
            nor.Normalize();
            normals[3*n+0] = nor.x;
            normals[3*n+1] = nor.y;
            normals[3*n+2] = nor.z;
            n++;
        }
    vertice[3*n+0] = 0;
    vertice[3*n+1] = -radius;
    vertice[3*n+2] = 0;
    normals[3*n+0] = 0;
    normals[3*n+1] = -1.0;
    normals[3*n+2] = 0;
    int w=sphereAxisSeg;
    int h=sphereHeightSeg;
    for (int i=0;i<w;i++) {
        face[i*3+0] = i+1;
        face[i*3+1] = 0;
        face[i*3+2] = (i+1)%w+1;
    }
    for (int j=0;j<h-2;j++)
        for (int i=0;i<w;i++) {
            face[(i+j*w)*6+w*3+0] = i+j*w+1;
            face[(i+j*w)*6+w*3+1] = (i+1)%w+j*w+1;
            face[(i+j*w)*6+w*3+2] = i+(j+1)*w+1;
            face[(i+j*w)*6+w*3+3] = (i+1)%w+j*w+1;
            face[(i+j*w)*6+w*3+4] = (i+1)%w+(j+1)*w+1;
            face[(i+j*w)*6+w*3+5] = i+(j+1)*w+1;
        }
    for (int i=0;i<w;i++) {
        face[w*(h-2)*6+w*3+i*3+0] = i+(h-2)*w+1;
        face[w*(h-2)*6+w*3+i*3+1] = (i+1)%w+(h-2)*w+1;
        face[w*(h-2)*6+w*3+i*3+2] = vertice_num-1;
    }
    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildCylinder(double x1, double y1, double z1, double radius1, double x2, double y2, double z2, double radius2)
{   
    type = CELL_CYLINDER;
    is_surface_constructed = TRIANGLES;
    if (vertice_num > 0) // to avoid double delete
    {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }

    // hard-coded parameters for now
    int numberOfSides = 50; // original 256

    cylinderRadius = radius1;

    double xDist = x2 - x1;
    double yDist = y2 - y1;
    double zDist = z2 - z1;
    cylinderHeight = sqrt(pow(xDist, 2) + pow(yDist, 2) + pow(zDist, 2));

    //MARK
    double dphi;
    cylinderAxisSeg = numberOfSides;
    dphi = 360.0/cylinderAxisSeg;       //aka there are cylinderAxisSeg number of wedges with angle dphi each
    cylinderHeightSeg = 1;

    vertice_num = (int)((cylinderHeightSeg+3)*360/dphi+2);
    face_num = (int)((360/dphi)*(cylinderHeightSeg))*2+(int)(360/dphi)*2;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];
    double theta;

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    int n = 0;  // index of vertices

    // center point of the top circle face
    vertice[3*n+0] = x1;
    vertice[3*n+1] = y1;
    vertice[3*n+2] = z1;
    normals[3*n+0] = 0;
    normals[3*n+1] = 1.0;
    normals[3*n+2] = 0;

    ++n; // increment 1 for center point of the top circle face

    // figure out basis vectors
    // e2 is parallel to the central axis
    Vector3 e2(x2 - x1, y2 - y1, z2 - z1);
    double e2norm = sqrt(pow(e2.x, 2) + pow(e2.y, 2) + pow(e2.z, 2));
    e2 = e2/(e2norm + DBL_EPSILON);

    // pick something -- anything that is perpendicular
    // pick something that isn't co-linear with e2
    Vector3 u1(-1 + 2*rand(), -1 + 2*rand(), -1 + 2*rand());

    // subtract off parallel part to get something fully perpendicular
    double dotProduct = u1.x * e2.x + u1.y * e2.y + u1.z * e2.z;
    Vector3 e1(u1.x - dotProduct*e2.x, u1.y - dotProduct*e2.y, u1.z - dotProduct*e2.z);
    // normalize it
    double e1norm = sqrt(pow(e1.x, 2) + pow(e1.y, 2) + pow(e1.z, 2));
    e1 = e1/e1norm;

    // pick a last basis vector perpendicular to both
    Vector3 e3;
    CellMath::crossProduct(e1.x, e1.y, e1.z, e2.x, e2.y, e2.z, e3.x, e3.y, e3.z);

    // points along the rim of the top circle face
    for (theta = 0; theta < 360 - CELL_EPS; theta += dphi)
    {
        vertice[3 * n] = x1 + cylinderRadius * (cosf(theta*DTOR) * e1.x + sinf(theta*DTOR) * e3.x);
        vertice[3 * n + 1] = y1 + cylinderRadius * (cosf(theta*DTOR) * e1.y + sinf(theta*DTOR) * e3.y);
        vertice[3 * n + 2] = z1 + cylinderRadius * (cosf(theta*DTOR) * e1.z + sinf(theta*DTOR) * e3.z);
        Vector3 norm(vertice[3 * n], vertice[3 * n + 1], vertice[3 * n + 2]);
        norm.Normalize();
        normals[3*n] = norm.x;
        normals[3*n + 1] = norm.y;
        normals[3*n + 2] = norm.z;
        ++n; // increment once for each point along the rim of the top circle face
    }

    // points for the top of the sides of the cylinder
    // note: this generates a cylinder for cylinderHeightSeg size of 1, aka not really any segmentation
    for (theta = 0; theta < 360 - CELL_EPS; theta += dphi)
    {
        vertice[3 * n] = x1 + cylinderRadius * (cosf(theta*DTOR) * e1.x + sinf(theta*DTOR) * e3.x);
        vertice[3 * n + 1] = y1 + cylinderRadius * (cosf(theta*DTOR) * e1.y + sinf(theta*DTOR) * e3.y);
        vertice[3 * n + 2] = z1 + cylinderRadius * (cosf(theta*DTOR) * e1.z + sinf(theta*DTOR) * e3.z);
        Vector3 norm(vertice[3 * n], vertice[3 * n + 1], vertice[3 * n + 2]);
        norm.Normalize();
        normals[3*n] = norm.x;
        normals[3*n + 1] = norm.y;
        normals[3*n + 2] = norm.z;
        ++n; // increment once for each point along the rim of the top circle face
    }

    // points for the bottom of the sides of the cylinder
    for (theta=0; theta < 360 - CELL_EPS; theta += dphi)
    {
        vertice[3 * n] = x2 + radius2*(cosf(theta*DTOR)*e1.x + sinf(theta*DTOR)*e3.x);
        vertice[3 * n + 1] = y2 + radius2*(cosf(theta*DTOR)*e1.y + sinf(theta*DTOR)*e3.y);
        vertice[3 * n + 2] = z2 + radius2*(cosf(theta*DTOR)*e1.z + sinf(theta*DTOR)*e3.z);
        Vector3 nor = Vector3(vertice[3*n], vertice[3*n+1], vertice[3*n+2]);
        nor.Normalize();
        normals[3*n+0] = nor.x;
        normals[3*n+1] = nor.y;
        normals[3*n+2] = nor.z;
        ++n;
    }

    //points along the rim of the bottom circle face
    for (theta=0; theta < 360 - CELL_EPS; theta += dphi)
    {
        vertice[3 * n] = x2 + radius2*(cosf(theta*DTOR)*e1.x + sinf(theta*DTOR)*e3.x);
        vertice[3 * n + 1] = y2 + radius2*(cosf(theta*DTOR)*e1.y + sinf(theta*DTOR)*e3.y);
        vertice[3 * n + 2] = z2 + radius2*(cosf(theta*DTOR)*e1.z + sinf(theta*DTOR)*e3.z);
        Vector3 nor = Vector3(vertice[3*n], vertice[3*n+1], vertice[3*n+2]);
        nor.Normalize();
        normals[3*n+0] = nor.x;
        normals[3*n+1] = nor.y;
        normals[3*n+2] = nor.z;
        ++n;
    }

    //center point for the bottom circle face
    vertice[3*n+0] = x2;
    vertice[3*n+1] = y2;
    vertice[3*n+2] = z2;
    normals[3*n+0] = 0;
    normals[3*n+1] = -1.0;
    normals[3*n+2] = 0;

    int w=360/dphi;
    for (int i=0;i<w;i++)
    {
        face[i*3+0] = i+1;
        face[i*3+1] = 0;
        face[i*3+2] = (i+1)%w+1;
    }
    for (int j = 0;j < cylinderHeightSeg; ++j) {
        for (int i=0;i<w;i++) {
            face[(i+j*w)*6+w*3+0] = i+j*w+1+w;
            face[(i+j*w)*6+w*3+1] = (i+1)%w+j*w+1+w;
            face[(i+j*w)*6+w*3+2] = i+(j+1)*w+1+w;
            face[(i+j*w)*6+w*3+3] = (i+1)%w+j*w+1+w;
            face[(i+j*w)*6+w*3+4] = (i+1)%w+(j+1)*w+1+w;
            face[(i+j*w)*6+w*3+5] = i+(j+1)*w+1+w;
        }
    }
    for (int i=0;i<w;i++)
    {
        face[w*(cylinderHeightSeg)*6+w*3+i*3+0] = i+cylinderHeightSeg*w+1+2*w;
        face[w*(cylinderHeightSeg)*6+w*3+i*3+1] = (i+1)%w+cylinderHeightSeg*w+1+2*w;
        face[w*(cylinderHeightSeg)*6+w*3+i*3+2] = vertice_num-1;
    }

    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildCylinder(double radius, double height, int wseg, int hseg)
{
    type = CELL_CYLINDER;
    is_surface_constructed = TRIANGLES;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }

    if (radius == 0.0)
    {
        radius = 5.0;
    }
    cylinderRadius = radius;

    if (height == 0.0)
    {
        height = 5.0;
    }
    cylinderHeight = height;

    double dphi;
    int hden;
    if (wseg==0)
    {
        dphi = 15;
        cylinderAxisSeg = 24;
    }
    else
    {
        dphi = 360.0/wseg;
        cylinderAxisSeg = wseg;
    }
    if (hseg==0)
    {
        hden = 1;
    }
    else
    {
        hden = hseg;
    }
    cylinderHeightSeg = hden;

    vertice_num = (int)((hden+3)*360/dphi+2);
    face_num = (int)((360/dphi)*(hden))*2+(int)(360/dphi)*2;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];
    double theta;

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    int n = 0;    //center point of the top circle face
    vertice[3*n+0] = 0;
    //vertice[3*n+1] = height/2;
    vertice[3*n+1] = height;
    vertice[3*n+2] = 0;
    normals[3*n+0] = 0;
    normals[3*n+1] = 1.0;
    normals[3*n+2] = 0;
    //points going along the rim of the top circle face
    n++;
    for (theta=0; theta<360-CELL_EPS; theta += dphi) {
        vertice[3*n+0] = radius*cosf(theta*DTOR);
        //vertice[3*n+1] = height/2;
        vertice[3*n+1] = height;
        vertice[3*n+2] = radius*sinf(DTOR*theta);
        normals[3*n+0] = 0;
        normals[3*n+1] = 1.0;
        normals[3*n+2] = 0;
        n++;
    }
    //points for the sides of the cylinder
    for (int i=0; i<=hden; i++) {
        for (theta=0; theta<360-CELL_EPS; theta += dphi) {
            vertice[3*n+0] = radius*cosf(theta*DTOR);
            vertice[3*n+1] = height - height/hden*i;
            //vertice[3*n+1] = height/2-height/hden*i;
            vertice[3*n+2] = radius*sinf(DTOR*theta);
            Vector3 nor = Vector3(vertice[3*n+0], 0.0, vertice[3*n+2]);
            nor.Normalize();
            normals[3*n+0] = nor.x;
            normals[3*n+1] = nor.y;
            normals[3*n+2] = nor.z;
            n++;
        }
    }
    //points for the bottom circle face
    for (theta=0; theta<360-CELL_EPS; theta += dphi) {
        vertice[3*n+0] = radius*cosf(theta*DTOR);
        //vertice[3*n+1] = -height/2;
        vertice[3*n+1] = 0;
        vertice[3*n+2] = radius*sinf(DTOR*theta);
        normals[3*n+0] = 0;
        normals[3*n+1] = -1.0;
        normals[3*n+2] = 0;
        n++;
    }
    //center point for the bottom circle face
    vertice[3*n+0] = 0;
    //vertice[3*n+1] = -height/2;
    vertice[3*n+1] = 0;
    vertice[3*n+2] = 0;
    normals[3*n+0] = 0;
    normals[3*n+1] = -1.0;
    normals[3*n+2] = 0;
    int w=360/dphi;
    for (int i=0;i<w;i++) {
        face[i*3+0] = i+1;
        face[i*3+1] = 0;
        face[i*3+2] = (i+1)%w+1;
    }
    for (int j=0;j<hden;j++)
        for (int i=0;i<w;i++) {
            face[(i+j*w)*6+w*3+0] = i+j*w+1+w;
            face[(i+j*w)*6+w*3+1] = (i+1)%w+j*w+1+w;
            face[(i+j*w)*6+w*3+2] = i+(j+1)*w+1+w;
            face[(i+j*w)*6+w*3+3] = (i+1)%w+j*w+1+w;
            face[(i+j*w)*6+w*3+4] = (i+1)%w+(j+1)*w+1+w;
            face[(i+j*w)*6+w*3+5] = i+(j+1)*w+1+w;
        }
    for (int i=0;i<w;i++) {
        face[w*(hden)*6+w*3+i*3+0] = i+hden*w+1+2*w;
        face[w*(hden)*6+w*3+i*3+1] = (i+1)%w+hden*w+1+2*w;
        face[w*(hden)*6+w*3+i*3+2] = vertice_num-1;
    }
    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildOpenCylinder(double radius, double height, int wseg, int hseg)
{
    type = CELL_OPENCYLINDER;
    is_surface_constructed = TRIANGLES;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }
    if (radius==0.0)
        radius = 5.0;
    cylinderRadius = radius;
    if (height==0.0)
        height = 5.0;
    cylinderHeight = height;

    double dphi;
    int hden;
    if (wseg==0) {
        dphi = 15;
        cylinderAxisSeg = 24;
    }
    else {
        dphi = 360.0/wseg;
        cylinderAxisSeg = wseg;
    }
    if (hseg==0)
        hden = 1;
    else
        hden = hseg;
    cylinderHeightSeg = hden;

    vertice_num = (int)((hden+1)*360/dphi);
    face_num = (int)((360/dphi)*(hden))*2;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];
    double theta;

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    int n = 0;
    for (int i=0; i<=hden; i++)
        for (theta=0; theta<360-CELL_EPS; theta += dphi) {
            vertice[3*n+0] = radius*cosf(theta*DTOR);
            vertice[3*n+1] = height/2-height/hden*i;
            vertice[3*n+2] = radius*sinf(DTOR*theta);
            Vector3 nor = Vector3(vertice[3*n+0], 0.0, vertice[3*n+2]);
            nor.Normalize();
            normals[3*n+0] = nor.x;
            normals[3*n+1] = nor.y;
            normals[3*n+2] = nor.z;
            n++;
        }

    int w=360/dphi;
    for (int j=0;j<hden;j++)
        for (int i=0;i<w;i++) {
            face[(i+j*w)*6+0] = i+j*w;
            face[(i+j*w)*6+1] = (i+1)%w+j*w;
            face[(i+j*w)*6+2] = i+(j+1)*w;
            face[(i+j*w)*6+3] = (i+1)%w+j*w;
            face[(i+j*w)*6+4] = (i+1)%w+(j+1)*w;
            face[(i+j*w)*6+5] = i+(j+1)*w;
        }

    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildPlane(double xlength, double zlength, int xseg, int zseg)
{
    type = CELL_PLANE;
    is_surface_constructed = TRIANGLES;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }

    //if width is 0, set it equal to 5.0
    if (xlength==0.0) {
        xlength = 5.0;
    }
    planeWidth = xlength;
    //if height is 0, set it equal to 5.0
    if (zlength==0.0) {
        zlength = 5.0;
    }
    planeLength = zlength;
    //if width segmentation is 0, set it equal to 1
    int xden, zden;
    if (xseg==0) {
        xden=1;
    }
    else {
        xden=xseg;
    }
    planeWidthSeg = xden;
    //if length segmentation is 0, set it equal to 1
    if (zseg==0) {
        zden=1;
    }
    else {
        zden=zseg;
    }
    planeLengthSeg = zden;
    //set number of vertices
    vertice_num = (xden+1)*(zden+1)*2;
    //set number of faces (not faces in the polygon, but segment faces, aka triangles)
    face_num = xden*zden*2;
    //allocate arrays to hold the vertices, normals, & faces
    //give extra space to vertice & normals
    vertice = new GLdouble[vertice_num*3];    
    normals = new GLdouble[vertice_num*3];
    //give exact space for face
    face = new int[face_num*3];

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    int n=0;
    int fn=0;
    for (int i=0; i<xden; i++)
        for (int j=0; j<zden; j++) {
            face[6*fn+0] = i+j*(xden+1);
            face[6*fn+1] = i+1+j*(xden+1);
            face[6*fn+2] = i+(j+1)*(xden+1);
            face[6*fn+3] = i+1+j*(xden+1);
            face[6*fn+4] = i+1+(j+1)*(xden+1);
            face[6*fn+5] = i+(j+1)*(xden+1);
            fn++;
        }
    for (int i=0;i<=xden;i++) {
        for (int j=0;j<=zden;j++)
        {
            /*
            //for plane in the xz-plane
            vertice[3*n+0] = xlength/xden*i-xlength/2;
            vertice[3*n+1] = 1.0;
            vertice[3*n+2] = zlength/zden*j-zlength/2;
            normals[3*n+0] = 0.0;
            normals[3*n+1] = 1.0;
            normals[3*n+2] = 0;
            */

            //for a plane in the yz-plane
            vertice[3*n+0] = 0.0;
            vertice[3*n+1] = xlength/xden*i-xlength/2;
            vertice[3*n+2] = zlength/zden*j-zlength/2;
            normals[3*n+0] = 1.0;
            normals[3*n+1] = 0;
            normals[3*n+2] = 0;
            n++;
        }
    }

    //reset deletePoint, selectedPoint, and deleteFace
    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildCube(double xlength, double ylength, double zlength, int xseg, int yseg, int zseg)
{
    type = CELL_CUBE;
    is_surface_constructed = TRIANGLES;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }
    if (xlength==0.0)
        xlength = 5.0;
    cubeWidth = xlength;
    if (ylength==0.0)
        ylength = 5.0;
    cubeHeight = ylength;
    if (zlength==0.0)
        zlength = 5.0;
    cubeDepth = zlength;
    int xden, yden, zden;
    if (xseg==0)
        xden=1;
    else
        xden=xseg;
    cubeWidthSeg = xden;
    if (yseg==0)
        yden=1;
    else
        yden=yseg;
    cubeHeightSeg = yden;
    if (zseg==0)
        zden=1;
    else
        zden=zseg;
    cubeDepthSeg = zden;

    vertice_num = (xden+1)*(yden+1)*2+(xden+1)*(zden+1)*2+(zden+1)*(yden+1)*2;
    face_num = (xden*yden*2+xden*zden*2+zden*yden*2)*2;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    int n=0;
    int fn=0;
    for (int i=0; i<xden; i++)
        for (int j=0; j<zden; j++) {
            face[6*fn+0] = j+i*(zden+1);
            face[6*fn+1] = j+1+i*(zden+1);
            face[6*fn+2] = j+1+(i+1)*(zden+1);
            face[6*fn+3] = j+i*(zden+1);
            face[6*fn+4] = j+1+(i+1)*(zden+1);
            face[6*fn+5] = j+(i+1)*(zden+1);
            fn++;
        }
    for (int i=0;i<=xden;i++)
        for (int j=0;j<=zden;j++) {
            vertice[3*n+0] = xlength/xden*i-xlength/2;
            vertice[3*n+1] = ylength/2;
            vertice[3*n+2] = zlength/zden*j-zlength/2;
            normals[3*n+0] = 0;
            normals[3*n+1] = 1.0;
            normals[3*n+2] = 0;
            n++;
        }
    for (int i=0; i<xden; i++)
        for (int j=0; j<zden; j++) {
            face[6*fn+0] = j+i*(zden+1)+n;
            face[6*fn+1] = j+1+(i+1)*(zden+1)+n;
            face[6*fn+2] = j+1+i*(zden+1)+n;
            face[6*fn+3] = j+i*(zden+1)+n;
            face[6*fn+4] = j+(i+1)*(zden+1)+n;
            face[6*fn+5] = j+1+(i+1)*(zden+1)+n;
            fn++;
        }
    for (int i=0;i<=xden;i++)
        for (int j=0;j<=zden;j++) {
            vertice[3*n+0] = xlength/xden*i-xlength/2;
            vertice[3*n+1] = -ylength/2;
            vertice[3*n+2] = zlength/zden*j-zlength/2;
            normals[3*n+0] = 0;
            normals[3*n+1] = -1.0;
            normals[3*n+2] = 0;
            n++;
        }
    for (int i=0; i<xden; i++)
        for (int j=0; j<yden; j++) {
            face[6*fn+0] = j+i*(yden+1)+n;
            face[6*fn+1] = j+1+(i+1)*(yden+1)+n;
            face[6*fn+2] = j+1+i*(yden+1)+n;
            face[6*fn+3] = j+i*(yden+1)+n;
            face[6*fn+4] = j+(i+1)*(yden+1)+n;
            face[6*fn+5] = j+1+(i+1)*(yden+1)+n;
            fn++;
        }
    for (int i=0;i<=xden;i++)
        for (int j=0;j<=yden;j++) {
            vertice[3*n+0] = xlength/xden*i-xlength/2;
            vertice[3*n+1] = ylength/yden*j-ylength/2;
            vertice[3*n+2] = zlength/2;
            normals[3*n+0] = 0.0;
            normals[3*n+1] = 0.0;
            normals[3*n+2] = 1.0;
            n++;
        }
    for (int i=0; i<xden; i++)
        for (int j=0; j<yden; j++) {
            face[6*fn+0] = j+i*(yden+1)+n;
            face[6*fn+1] = j+1+i*(yden+1)+n;
            face[6*fn+2] = j+1+(i+1)*(yden+1)+n;
            face[6*fn+3] = j+i*(yden+1)+n;
            face[6*fn+4] = j+1+(i+1)*(yden+1)+n;
            face[6*fn+5] = j+(i+1)*(yden+1)+n;
            fn++;
        }
    for (int i=0;i<=xden;i++)
        for (int j=0;j<=yden;j++) {
            vertice[3*n+0] = xlength/xden*i-xlength/2;
            vertice[3*n+1] = ylength/yden*j-ylength/2;
            vertice[3*n+2] = -zlength/2;
            normals[3*n+0] = 0.0;
            normals[3*n+1] = 0.0;
            normals[3*n+2] = -1.0;
            n++;
        }
    for (int i=0; i<zden; i++)
        for (int j=0; j<yden; j++) {
            face[6*fn+0] = j+i*(yden+1)+n;
            face[6*fn+1] = j+1+i*(yden+1)+n;
            face[6*fn+2] = j+1+(i+1)*(yden+1)+n;
            face[6*fn+3] = j+i*(yden+1)+n;
            face[6*fn+4] = j+1+(i+1)*(yden+1)+n;
            face[6*fn+5] = j+(i+1)*(yden+1)+n;
            fn++;
        }
    for (int i=0;i<=zden;i++)
        for (int j=0;j<=yden;j++) {
            vertice[3*n+0] = xlength/2;
            vertice[3*n+1] = ylength/yden*j-ylength/2;
            vertice[3*n+2] = zlength/zden*i-zlength/2;
            normals[3*n+0] = 1.0;
            normals[3*n+1] = 0.0;
            normals[3*n+2] = 0.0;
            n++;
        }
    for (int i=0; i<zden; i++)
        for (int j=0; j<yden; j++) {
            face[6*fn+0] = j+i*(yden+1)+n;
            face[6*fn+1] = j+1+(i+1)*(yden+1)+n;
            face[6*fn+2] = j+1+i*(yden+1)+n;
            face[6*fn+3] = j+i*(yden+1)+n;
            face[6*fn+4] = j+(i+1)*(yden+1)+n;
            face[6*fn+5] = j+1+(i+1)*(yden+1)+n;
            fn++;
        }
    for (int i=0;i<=zden;i++)
        for (int j=0;j<=yden;j++) {
            vertice[3*n+0] = -xlength/2;
            vertice[3*n+1] = ylength/yden*j-ylength/2;
            vertice[3*n+2] = zlength/zden*i-zlength/2;
            normals[3*n+0] = -1.0;
            normals[3*n+1] = 0.0;
            normals[3*n+2] = 0.0;
            n++;
        }

    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildCone(double height, double radius, int wseg)
{
    is_surface_constructed = TRIANGLES;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }
    if (radius==0.0)
        radius = 5.0;
    if (height==0.0)
        height = 10.0;
    double dphi;
    if (wseg==0)
        dphi = 5;
    else
        dphi = 360.0/wseg;
    vertice_num = (int)(360/dphi*2+2);
    face_num = (int)(360/dphi)*2;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    int n = 0;
    vertice[3*n+0] = 0;
    vertice[3*n+1] = height;
    vertice[3*n+2] = 0;
    normals[3*n+0] = 0;
    normals[3*n+1] = 1.0;
    normals[3*n+2] = 0;
    n++;
    for (int theta=0; theta<360; theta += (int)dphi) {
        vertice[3*n+0] = radius*cosf(theta*DTOR);
        vertice[3*n+1] = 0.0;
        vertice[3*n+2] = radius*sinf(DTOR*theta);
        Vector3 nor = Vector3(vertice[3*n+0], 0.0, vertice[3*n+2]);
        nor.Normalize();
        normals[3*n+0] = nor.x;
        normals[3*n+1] = nor.y;
        normals[3*n+2] = nor.z;
        n++;
    }
    for (int theta=0; theta<360; theta += (int)dphi) {
        vertice[3*n+0] = radius*cosf(theta*DTOR);
        vertice[3*n+1] = 0.0;
        vertice[3*n+2] = radius*sinf(DTOR*theta);
        normals[3*n+0] = 0;
        normals[3*n+1] = -1.0;
        normals[3*n+2] = 0;
        n++;
    }
    vertice[3*n+0] = 0;
    vertice[3*n+1] = 0;
    vertice[3*n+2] = 0;
    normals[3*n+0] = 0;
    normals[3*n+1] = -1.0;
    normals[3*n+2] = 0;
    int w=360/dphi;
    for (int i=0;i<w;i++) {
        face[i*3+0] = i+1;
        face[i*3+1] = 0;
        face[i*3+2] = (i+1)%w+1;
    }
    for (int i=w;i<w*2;i++) {
        face[i*3+0] = i-w+1+w;
        face[i*3+1] = (i-w+1)%w+1+w;
        face[i*3+2] = w+1+w;
    }
    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildLine(double ylength)
{
    type = CELL_LINE;
    is_surface_constructed = LINE_SEGMENT;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }

    if (ylength==0.0) {
        ylength = 5.0;
    }

    vertice_num = 2;
    face_num = 0;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    vertice[0] = 0.0;
    vertice[1] = -ylength/2;
    vertice[2] = 0.0;
    vertice[3] = 0.0;
    vertice[4] = ylength/2;
    vertice[5] = 0.0;

    //reset deletePoint, selectedPoint, and deleteFace
    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::buildObject(int v_num, int f_num, double *ver, double *nor, int *f)
{
    type = CELL_MESH;
    is_surface_constructed = TRIANGLES;
    if (vertice_num>0) {
        delete [] vertice;
        delete [] normals;
        delete [] face;
    }
    vertice_num = v_num;
    face_num = f_num;
    vertice = new GLdouble[vertice_num*3];
    normals = new GLdouble[vertice_num*3];
    face = new int[face_num*3];

    transformedVertices = new GLdouble[vertice_num*3];
    transformedNormals = new GLdouble[vertice_num*3];
    transformedFaces = new int[face_num*3];

    for (int i=0;i<vertice_num*3;i++) {
        vertice[i] = ver[i];
        normals[i] = nor[i];
    }
    for (int i=0;i<face_num*3;i++) {
        face[i] = f[i];
    }
    deletePoint = new bool[vertice_num];
    selectedPoint = new bool[vertice_num];
    resetSelectedPoint();
    for (int i=0;i<vertice_num;i++)
        deletePoint[i]=false;
    deleteFace = new bool[face_num];
    for (int i=0;i<face_num;i++)
        deleteFace[i]=false;
    selectedFace = new bool[face_num];
    resetSelectedFace();
}

void CellObject::toCSG()
{
    if (isCSG == false) {
        csgObject = new CSG();
        transformVertex();
        vector<MyPolygon> pa;
        for (int i=0;i<face_num;i++) {
            Vertex3 v0 = Vertex3(this->getVertice(this->getFace(i*3+0)), this->getNormal(this->getFace(i*3+0)));
            Vertex3 v1 = Vertex3(this->getVertice(this->getFace(i*3+1)), this->getNormal(this->getFace(i*3+1)));
            Vertex3 v2 = Vertex3(this->getVertice(this->getFace(i*3+2)), this->getNormal(this->getFace(i*3+2)));
            vector<Vertex3> v;
            v.push_back(v0);
            v.push_back(v1);
            v.push_back(v2);
            MyPolygon p = MyPolygon(v, 0.0);
            pa.push_back(p);
        }
        csgObject->polygons = pa;
    } else {
        transformCSG();
    }
}

void CellObject::setCSG(CSG *csg)
{
    csgObject = new CSG();
    csgObject->polygons = csg->polygons;
    isCSG = true;
}

void CellObject::toObject()
{
    int f_num = 0;
    int v_num = 0;
    for (size_t i=0;i<csgObject->polygons.size();i++) {
        f_num += (csgObject->polygons[i]).vertices.size()-2;
        v_num += (csgObject->polygons[i]).vertices.size();
    }
    int *face = new int[f_num*3];
    double *ver = new double[v_num*3];
    double *nor = new double[v_num*3];
    int fn = 0;
    int vn = 0;
    for (size_t i=0;i<csgObject->polygons.size();i++) {
        int tmpvn = (csgObject->polygons[i]).vertices.size();
        for (int k=2;k<tmpvn;k++) {
            face[fn*3+0] = vn;
            face[fn*3+1] = vn+k-1;
            face[fn*3+2] = vn+k;
            fn++;
        }
        for (int k=0;k<tmpvn;k++) {
            ver[vn*3+0] = (csgObject->polygons[i]).vertices[k].pos.x;
            ver[vn*3+1] = (csgObject->polygons[i]).vertices[k].pos.y;
            ver[vn*3+2] = (csgObject->polygons[i]).vertices[k].pos.z;
            nor[vn*3+0] = (csgObject->polygons[i]).vertices[k].normal.x;
            nor[vn*3+1] = (csgObject->polygons[i]).vertices[k].normal.y;
            nor[vn*3+2] = (csgObject->polygons[i]).vertices[k].normal.z;
            vn++;
        }
    }
    buildObject(v_num, f_num, ver, nor, face);
    delete [] face;
    delete [] ver;
    delete [] nor;
}

void CellObject::addRotateX(double x)
{
    double matX[4][4], result[4][4];
    CellMath::rotateMatX(x, matX);
    CellMath::matMutiply(matX, newRotateMatrix, result);

    for(int i = 0; i < 4; i++) {
        for(int j = 0; j < 4; j++) {
            newRotateMatrix[i][j] = result[i][j];
        }
    }
}

void CellObject::addRotateY(double y)
{
    double matY[4][4], result[4][4];
    CellMath::rotateMatY(y, matY);
    CellMath::matMutiply(matY, newRotateMatrix, result);

    for(int i = 0; i < 4; i++) {
        for(int j = 0; j < 4; j++) {
            newRotateMatrix[i][j] = result[i][j];
        }
    }
}

void CellObject::addRotateZ(double z)
{
    double matZ[4][4], result[4][4];
    CellMath::rotateMatZ(z, matZ);
    CellMath::matMutiply(matZ, newRotateMatrix, result);

    for(int i = 0; i < 4; i++) {
        for(int j = 0; j < 4; j++) {
            newRotateMatrix[i][j] = result[i][j];
        }
    }
}

void CellObject::addRotateCtrl(double a, double x, double y, double z)
{
//    double rotMat[4][4], result[4][4];
//    CellMath::rotateMatFromAxis(a, x, y, z, rotMat);
//    CellMath::matMutiply(rotMat, newRotateMatrix, result);

//    for(int i = 0; i < 4; i++) {
//        for(int j = 0; j < 4; j++) {
//            newRotateMatrix[i][j] = result[i][j];
//        }
//    }
    rotateCtrl += a;
    rotateAxisX = x;
    rotateAxisY = y;
    rotateAxisZ = z;
}

void CellObject::addTranslateCtrlPt1X(double x)
{
    if (type == CELL_LINE)
    {
        vertice[0] += x;
    }
    else
    {
        translateX += x;
    }
}

void CellObject::addTranslateCtrlPt1Y(double y)
{
    if (type == CELL_LINE)
    {
        vertice[1] += y;
    }
    else
    {
        translateY += y;
    }
}

void CellObject::addTranslateCtrlPt1Z(double z)
{
    if (type == CELL_LINE)
    {
        vertice[2] += z;
    }
    else
    {
        translateZ += z;
    }
}

void CellObject::addTranslateCtrlPt2X(double x)
{
    if (type == CELL_LINE)
    {
        vertice[3] += x;
    }
    else
    {
        translateX += x;
    }
}

void CellObject::addTranslateCtrlPt2Y(double y)
{
    if (type == CELL_LINE)
    {
        vertice[4] += y;
    }
    else
    {
        translateY += y;
    }
}

void CellObject::addTranslateCtrlPt2Z(double z)
{
    if (type == CELL_LINE)
    {
        vertice[5] += z;
    }
    else
    {
        translateZ += z;
    }
}

void CellObject::addTranslateX(double x)
{
    translateX += x;
}

void CellObject::addTranslateY(double y)
{
    translateY += y;
}

void CellObject::addTranslateZ(double z)
{
    translateZ += z;
}
